#!/bin/bash
#从 kingdom 导出层位然后 line drawing,由于 horizon 需要新建很多层位,因此用断层表示层位,导出的都是断层格式(GeoQuest-CARD IMAGE 7),注意左对齐
function preset()
{
    #=================================
    it=0
    line=04
    fig_name=${line}_normalfault_drawing1
    fig_fmt=png
    dmig=../../../sgy/grd/${line}_migkd_it${it}.grd
    #=================================
    #gmt set MAP_FRAME_TYPE=plain
    gmt makecpt -Cwgb.cpt -T-10000/10000 > wgb_25.cpt
    MC32=240021.54
    MC11=256856.42
    MC15=270377.61
    MC16A=287319.20
    MC30=295054.90
    MC13=318341.85
    
}

preset

gmt begin $fig_name $fig_fmt E300 
# ============================================================ (a) seismic profile===========================================================================================
origin_x_min=9200 #set min_cmp num 
origin_x_max=10500 #set max_cmp num 
new_x_min=`echo "scale=2; $origin_x_min * 6.25 / 1000 " | bc`
new_x_max=`echo "scale=2; $origin_x_max * 6.25 / 1000 " | bc`
origin_d_min=4.741
origin_d_max=9.884
xylimits_cmp=`echo $origin_x_min/$origin_x_max/$origin_d_min/$origin_d_max`
trench_cmp=6900
trench_km=`echo "scale=2; $trench_cmp * 6.25 / 1000 " | bc`
left_trench=`echo "scale=2; $new_x_min - $trench_km " | bc`
right_trench=`echo "scale=2; $new_x_max - $trench_km " | bc`
xylimits_km=`echo $left_trench/$right_trench/$origin_d_min/$origin_d_max` #6900*6.25/1000
height=`echo "scale=2; ($origin_d_max - $origin_d_min)/3" | bc`
vert_exa=1
width=`echo "scale=2; $height * (${origin_x_max} - ${origin_x_min}) * 6.25 *0.001  / ((${origin_d_max} - ${origin_d_min})* ${vert_exa})" | bc`
y_shift=`echo "scale=2; $width+0.3 " | bc`
x_shift=`echo "scale=2; $height+0.8 " | bc`
for v in ${line}  #读出起始终止坐标
do
x1=$(awk '{if($1=='$origin_x_min') print $2}' ../../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
y1=$(awk '{if($1=='$origin_x_min') print $3}' ../../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
x2=$(awk '{if($1=='$origin_x_max') print $2}' ../../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
y2=$(awk '{if($1=='$origin_x_max') print $3}' ../../../gmt/P1_basemap_shade/MGL1610MC${v}_geometry_grd_rot_crooked.axon_MC${v}cmp_xy_m.txt) #utm
done
gmt grdmath $dmig 0.016 MUL = tmp.grd
gmt grdimage tmp.grd -Cwgb_25.cpt -R$xylimits_cmp -JX${width}c/-${height}c -B+n -M
gmt basemap -R$x1/$x2/$origin_d_min/$origin_d_max -JX${width}c/-${height}c -Byf0.5a1.0+l"Depth (km)" -BWesn  \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p --MAP_FRAME_PEN=0.5p --MAP_ANNOT_OFFSET=1p
gmt basemap -R$xylimits_km -JX${width}c/-${height}c -Bxf5a5+l"Distance to trench (km)" -BS \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p --MAP_FRAME_PEN=0.5p --MAP_ANNOT_OFFSET=1p
gmt basemap -R$xylimits_cmp -JX${width}c/-${height}c -Bxf100a2000+l"CMP number" -BN \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p --MAP_FRAME_PEN=0.5p --MAP_ANNOT_OFFSET=1p
gmt text -F+cBL+f5p,1 -Dj0.1c/0.1c <<EOF
MC${line}
EOF
gmt text -F+cTC+f5p,1 -Dj0.1c/0.1c <<EOF
A
EOF
gmt text -F+cBR+f5p,1 -Dj0.1c/0.1c <<EOF
VE=1
EOF
gmt text -F+cTL+f5p,1 -Dj0.1c/0.1c <<EOF
W
EOF
gmt text -F+cTR+f5p,1 -Dj0.1c/0.1c <<EOF
E
EOF
#<<COMMENT
# ============================================================ (b) line drawing===========================================================================================

for v in ${line}  #seafloor
do
gmt basemap -R$x1/$x2/$origin_d_min/$origin_d_max -JX${width}c/-${height}c -Byf0.5a1.0+l"Depth (km)" -BwEsn -X${y_shift} \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p --MAP_FRAME_PEN=0.5p --MAP_ANNOT_OFFSET=1p
awk '{print $1,$3/1000*-1}' ../../slice_seismic/seafloor_val_${v}.utm > tmp.txt
gmt plot tmp.txt -W0.1p,black  
done
echo ${MC32} `echo "scale=2; ${origin_d_min} " | bc` > tmp.txt
echo ${MC32} `echo "scale=2; ${origin_d_max} " | bc` >> tmp.txt
gmt text -F+f5p,1 <<EOF
${MC32} `echo "scale=2; ${origin_d_max}-0.5 " | bc` MC32
EOF
echo ${MC11} `echo "scale=2; ${origin_d_min} " | bc` > tmp.txt
echo ${MC11} `echo "scale=2; ${origin_d_max} " | bc` >> tmp.txt
gmt text -F+f5p,1 <<EOF
${MC11} `echo "scale=2; ${origin_d_max}-0.5 " | bc` MC11
EOF
echo ${MC15} `echo "scale=2; ${origin_d_min} " | bc` > tmp.txt
echo ${MC15} `echo "scale=2; ${origin_d_max} " | bc` >> tmp.txt
gmt text -F+f5p,1 <<EOF
${MC15} `echo "scale=2; ${origin_d_max}-0.5 " | bc` MC15
EOF
echo ${MC16A} `echo "scale=2; ${origin_d_min} " | bc` > tmp.txt
echo ${MC16A} `echo "scale=2; ${origin_d_max} " | bc` >> tmp.txt
gmt text -F+f5p,1 <<EOF
${MC16A} `echo "scale=2; ${origin_d_max}-0.5 " | bc` MC16A
EOF
echo ${MC30} `echo "scale=2; ${origin_d_min} " | bc` > tmp.txt
echo ${MC30} `echo "scale=2; ${origin_d_max} " | bc` >> tmp.txt
gmt text -F+f5p,1 <<EOF
${MC30} `echo "scale=2; ${origin_d_max}-0.5 " | bc` MC30
EOF
echo ${MC13} `echo "scale=2; ${origin_d_min} " | bc` > tmp.txt
echo ${MC13} `echo "scale=2; ${origin_d_max} " | bc` >> tmp.txt
gmt text -F+f5p,1 <<EOF
${MC13} `echo "scale=2; ${origin_d_max}-0.5 " | bc` MC13
EOF
awk '$3~/[[:digit:]+]/{print $1,$4/1000,$3}' ../../line_drawing/${line}_basement.dat > tmp0.txt #basement
k=$(awk 'END{print $3}' tmp0.txt) #the max num of the layers
i=$(awk 'NR==1{print $3}' tmp0.txt)
while [ ${i} -le ${k} ]; do
awk '$3 == '${i}'{print $1,$2}' tmp0.txt > tmp.txt
gmt plot tmp.txt -W0.1p,black
i=$(($i + 1))
done

awk '$3~/[[:digit:]+]/{print $1,$4/1000,$3}' ../../line_drawing/${line}_horizon1.dat > tmp0.txt #horizon
k=$(awk 'END{print $3}' tmp0.txt) #the max num of the layers
i=$(awk 'NR==1{print $3}' tmp0.txt)
while [ ${i} -le ${k} ]; do
awk '$3 == '${i}'{print $1,$2}' tmp0.txt > tmp.txt
gmt plot tmp.txt -W0.1p,magenta
i=$(($i + 1))
done

awk '$3~/[[:digit:]+]/{print $1,$4/1000,$3}' ../../line_drawing/${line}_horizon2.dat > tmp0.txt #horizon
k=$(awk 'END{print $3}' tmp0.txt) #the max num of the layers
i=$(awk 'NR==1{print $3}' tmp0.txt)
while [ ${i} -le ${k} ]; do
awk '$3 == '${i}'{print $1,$2}' tmp0.txt > tmp.txt
gmt plot tmp.txt -W0.1p,magenta
i=$(($i + 1))
done

awk '$3~/[[:digit:]+]/{print $1,$4/1000,$3}' ../../line_drawing/${line}_f.dat > tmp0.txt #fault
k=$(awk 'END{print $3}' tmp0.txt) #the max num of the layers
i=$(awk 'NR==1{print $3}' tmp0.txt)
while [ ${i} -le ${k} ]; do
awk '$3 == '${i}'{print $1,$2}' tmp0.txt > tmp.txt
gmt plot tmp.txt -W0.1p,red
i=$(($i + 1))
done

awk '$3~/[[:digit:]+]/{print $1,$4/1000,$3}' ../../line_drawing/${line}_rotat.dat > tmp0.txt #rotat
k=$(awk 'END{print $3}' tmp0.txt) #the max num of the layers
i=$(awk 'NR==1{print $3}' tmp0.txt)
while [ ${i} -le ${k} ]; do
awk '$3 == '${i}'{print $1,$2}' tmp0.txt > tmp.txt
gmt plot tmp.txt -W0.1p,grey
i=$(($i + 1))
done

gmt basemap -R$xylimits_km -JX${width}c/-${height}c -Bxf5a5+l"Distance to trench (km)" -BS \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p --MAP_FRAME_PEN=0.5p --MAP_ANNOT_OFFSET=1p
gmt basemap -R$xylimits_cmp -JX${width}c/-${height}c -Bxf100a2000+l"CMP number" -BN \
--MAP_TICK_LENGTH_PRIMARY=-0.03i --FONT_ANNOT_PRIMARY=5p --FONT_LABEL=5p --MAP_LABEL_OFFSET=1p --MAP_FRAME_PEN=0.5p --MAP_ANNOT_OFFSET=1p
gmt text -F+cBL+f5p,1 -Dj0.1c/0.1c <<EOF
MC${line}
EOF
gmt text -F+cTC+f5p,1 -Dj0.1c/0.1c <<EOF
A
EOF
gmt text -F+cBR+f5p,1 -Dj0.1c/0.1c <<EOF
VE=1
EOF
gmt text -F+cTL+f5p,1 -Dj0.1c/0.1c <<EOF
W
EOF
gmt text -F+cTR+f5p,1 -Dj0.1c/0.1c <<EOF
E
EOF
#COMMENT
gmt end show
rm tmp*
rm gmt*